Modeling the formation of in vitro filopodia 
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Filopodia are bundles ol actin filaments that extend out ahead of the leading edge of a crawling 
cell to probe its upcoming environment. In vitro experiments [D. Vignjevic et al, J. Cell Biol. 
160, 951 (2003)] have determined the minimal ingredients required for the formation of filopodia 
from the dendritic-like morphology of the leading edge. We model these experiments using kinetic 
aggregation equations for the density of growing bundle tips. In mean field, we determine the 
bundle size distribution to be broad for bundle sizes smaller than a characteristic bundle size above 
which the distribution decays exponentially. Two-dimensional simulations incorporating both 
bundling and cross-linking measure a bundle size distribution that agrees qualitatively with mean 
field. The simulations also demonstrate a nonmonotonicity in the radial extent of the dendritic 
region as a function of capping protein concentration, as was observed in experiments, due to the 
interplay between percolation and the ratcheting of growing filaments off a spherical obstacle. 
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I. INTRODUCTION 



Just as our skeletons give shape to our bodies, the cytoskeleton gives shape to the cell. However, unlike our 
skeletons, the "bones" , or filaments, of the cytoskeleton are dynamic — lengthening, shortening, branching, etc. — so 
that the cell can extend itself to crawl along surfaces to search for food. The physical process of cell crawling consists 
of (1) formation of a broad, thin protrusion of actin filaments along part of the cell periphery (a lamellipod), (2) the 
assembly of adhesion sites adhering the protrusion to the surface and disassembly of adhesion sites in the "rear" of 
the cell, and (3) retraction of the rear to catch up with the "front" part where the new growth has occurred. 

To date there are a number of quantitative models addressing each aspect of cell motility [l|-[3| . For example, nu- 
merous observations of lamellipodia have led researchers to conjecture that lamcllipodia are initiated by the branching 
agent Arp2/3, which nucleates new filaments from prexisting ones at an angle of about 70 degrees from the growing 
end of the polar actin filaments The resulting structure forms a directed, branched/dendritic-like network of 

filaments at the leading edge. Quantitative models implementing this mechanism of formation are able to account for 
observed cell speeds and various observed morphological properties (§-H3|. 

Here we quantitatively investigate a particular aspect of this emergent protrusion, namely the formation of bundles 
of actin filaments that extend out ahead of lamellipodia [i~3l ]. These bundles are known as filopodia. It has been 
proposed that filopodia serve as scouts for the cell [14j . Given this important role for filopodia in cell crawling, we ask 
the following questions: How do filopodia form? More specifically, how do they emerge from a branched structure? 
What is the typical number of actin filaments in a filopod? To build a more comprehensive framework of filopodia 
formation, how do oriented bundles of actin filaments emerge in the absence of lamcllipodia? 

To begin to address the formation of filopodia, in vitro experiments conducted by Vignjevic and collaborators [i~5| , 
along with experiments followed up by Haviv and collaborators [l6[ and Ideses and collaborators [TtJ , have determined 
the minimal ingredients required for filopodia formation. In the initial experiment [lH ]. a bead coated with Arp2/3 
activator was inserted into a solution containing actin monomers, linker/bundling proteins, capping protein, and 
Arp2/3. After some minutes, for a large enough concentration of Arp2/3, actin filaments emerge and form a mesh-like 
phase close to the bead. Under certain conditions, bundles of parallel filaments form further away from the bead — 
resulting in a star-like formation. The more recent in vitro experiments by Haviv and collaborators are conducted 
without a bead coated with Arp2/3 activator, but with activator in solution. They find similar star formation under 
certain conditions. 

Svitkina and collaborators fi~c| proposed two mechanisms for filopodia formation given the findings of the initial 
experiment. The first is that filopodia are the result of cone- like structures, termed lambda precursors, consisting 



§ E-mail address of corresponding author: jschwarz@physics.syr.edu 



2 



of filaments converging as they emerge from the dendritic network and continue on to grow and bundle. Similar 
structures have also been observed in in vivo B16F1 melanoma cells Svitkina and collaborators [l8[ have dubbed 
this mechanism of formation, the convergent elongation model. A second mechanism consists of longer bundles 
coalescing laterally via zippering to form even larger bundles. This latter mechanism does not preclude the former, 
however, the latter is probably less likely given the larger-scale rearrangements required to bring two pre-formed 
bundles, each with their respective base structure, together. 

The notion of filaments elongating and bundling is a rather natural one beggingg for quantitative modelling. 
We use the powerful framework of kinetic aggregation, initiated by the 1916 work of Smoluchowski studying the 
size distribution of "sticky" Brownian particles [20j . to build a model for filopodia formation [2l|. We construct 
an equation of motion for filament tips that grow and bundle upon colliding as dictated by simple geometry with 
some caveats regarding the details of the bundling. We use the equation of motion for filament tips to compute 
the distribution of bundle sizes, a measureable quantity. No ener geti cs is incorporated into the model explicitly, an 
assumption that is supported by work by Yang and collaborators [22[ , where careful simulations of two constrained- 
at-one-end semiflexible polymers fluctuating and coalescing show that energetics denominates beyond a time scale of 
10 seconds. In the star system, filaments lengthen and nucleate on a much faster time scale and so we expect the star 
formation to be kinetically driven. 

There have been several quantitative kinetic modeling papers addressing the formation of filopodia in the context of 
these experiments. For example, Mogilner and Rubinstein |23[ investigate the dynamics of lambda precursors colliding 
and maturing into full-blown filopodia. They find that density of filopodia at the leading edge is proportional to the 
square root of the initiation rate of lambda precursors. More recent work by Kraikivski and collaborators (2~H 
investigates the tip-tip bundling events versus tip-side bundling events to find that, contrary to expectation, both 
types of events produce a comparable number of collisions. However, the number of bundling events scales differently 
with the average bundle length for each respective event, as expected. 

The previous two theoretical works do not address a bundle size distribution. There exists a bundle size distribution 
calculation in a work, inspired by the initial experiments, treating the growing filaments as Potts-like spins on a 
lattice [25| . However, the kinetic Monte Carlo analysis yields a bundle size distribution with a prominent peak near 
a bundle size of three filaments (spins). This number is somewhat small compared with the qualitative estimates 
from the experiments. In the absence of Arp2/3, quantitative in vitro measurements of the bundle size distribution 
show that the average bundle size ranges from 10-30 depending on the ratio of actin to cross-linker concentration [26| . 
Ref. [26j | also proposes a model for the bundle size distribution where there is a rate for filaments leaving the bundle as 
well as contributing to the bundle. Since the filaments are entangled, extended objects, the rate for leaving a bundle 
seems somewhat small, however. 

In this work, we calculate mean field and two-dimensional bundle size distributions, as opposed to the total number 
of bundling events. In mean field, we find that tip-tip bundling alone is not sufficient to generate substantive bundle 
sizes, i.e. an average bundle size of at least 10. In two-dimensions, we find a threshold concentration of filaments, 
unbranched or branched, for substantive bundling. In mean field and in two-dimensions we find that the bundle size 
distribution is broad for bundle sizes below some characteristic bundle size and exponentially decays for bundle sizes 
larger than the characteristic bundle size. 

We also probe the change in morphology from the dendritic "phase" to bundled "phase" and find that to observe a 
highly oriented/bundled phase (1) Arp2/3 needs to be excluded from the periphery of the system and (2) the bundling 
proteins must reorient the filaments/bundles. We then investigate the effects of capping protein on the star formation 
and find a nonomonotonic dependence on the radial distribution of filament tips as a function capping protein in the 
presence of the bead. The bead is an obstacle which filaments growing into the bead can ratchet off of. The higher 
the concentration of capping protein, the less likely a percolating structure will form across the system such that 
filaments, and their cross-linked and bundled neighbors, can ratchet off the bead and expand the dendritic region. 

The organization of the paper is as follows. Section II constructs a general model based on the kinetic aggregation 
of bundle tips to calculate the bundle size distribution. Section III addresses the model within mean field. Section IV 
introduces and tests a kinetic simulation to go beyond mean field. Section V discusses the implications of our results. 



II. MODEL 



Given the in vitro experiments, let us assume that there are two processes in the lamcllipodia-filopodia transition: 
(1) branching by Arp2/3 and (2) bundling by cross-linking proteins. Branching allows for the nucleation of filaments, 
while bundling coalesces smaller bundles into larger ones. We will also investigate the effects of capping to more 
closely model the initial in vitro filopodia experiment where the concentration of capping protein was varied. Capping 
results in "death" of a growing filament tip. These three processes are represented diagrammatically in Figure 1. In 
the branching process, a new filament emerges from a prexisting one. For the bundling/coalescing process, we assume 
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that when two growing tips intersect with an angle less than 9 C , the two bundles realign and bundle with probability 
unity, otherwise, no bundle forms [22], [53, [H| . In principle, for intersection angles larger than 9 C , a cross- link can still 
form without bundle reorientation. This large-angle cross-link will be addressed in Section IV. There are other ways 
to bundle. In particular, the growing end of one bundle could collide into the side of another bundle. We call this 
tip-side bundling. See Fig. l(e)-(f). We do not consider zippering because it is a slower, less likely process, which 
presumably depends on the dynamics of the cross-linkers. We also do not consider depolymerization as it occurs on 
a much slower time scale than polymerization. 

Both tip-tip bundling and capping occurs at bundle/filament tips. We also assume that branching occurs at the 
growing bundle/filament end as well. The assumption of branching at the tip underestimates the effects of branching 
since side branching has been rather well-established experimentally (29[. However, there is also a known aging 
effect where Arp2/3 is found to bind preferentially to ATP-actin in filament s 13011 . In other words, branching is 
biased towards the growing end, which has also been observed experimentally [3ll433l |. These assumptions allow us 
to introduce the variable, pi(x, t), the density of growing tips in a bundle of size i at position x and time t. Tip-side 
bundling can be incorporated into the model by introducing an average growing length of bundles. We will address 
this extension in the next section. 

How does one describe the dynamics of the growing tips? Since actin is a scmiflexiblc polymer with a persistence 
lengthscale of about 17 microns, and the experiments indicate a star size on the order of 10 microns, the tip itself docs 
not diffuse, but rather the angle of the tip with respect to some fixed axis does. In other words, over the lengthscale 
of the experiment, the motion of the filament tips is essentially ballistic. 

Putting all of the above ingredients together yields an equation of motion for p\ , 

<9 t pi(x,i) = -w p Vpi(x, t) + £(x)pi (x, t) + b/o^(t)o-(x, t) }~] gtpjjx, t) - cpi(x, t) + 2cp 2 (x, t) - 2Ai,i(x,t)p? (x,t), (1) 

i=l 

where £(x) accounts for the angular diffusion of the tip and v p is the polymerization velocity. We have allowed for 
spatial variations in the density of Arp2/3 denoted by cr(x, t). The density of G-actin is denoted by p m (x,t), It is 
known that the nucleation of a branch at rate b requires a dimer of G-actin monomers hence the quadratic dependence 
on p m (x,i). The gi factor is to account for any geometric factors related to the dimensionality of the system and the 
size of the bundle. Note that the branching depends on the total concentration (or number of) bundles as opposed 
to the total concentration of filament tips. Finally, capping occurs at some uniform rate c and the bundling of two 
single filaments occurs at some rate Ai i i(x, t). 

Since we have introduced the concentrations of Arp2/3 and G-actin, we now construct their respective equations 
of motion. Arp2/3 diffuses and get taken up into branches. Therefore, the density of Arp2/3 is given by 

d t a(x,t) = D Arp2/3 V 2 a(x,t) -bp 2 m (t)a(x:,t)J29iPi(x,t), (2) 

i=l 

where .0^4^2/3 is the diffusion coefficient for Arp2/3. The concentration of G-actin is governed by 

d t p m {x,t) = D m V 2 p m {x,t) - k p ipi (x, t)p m (x, t) - 6p^(x, t)a(x, t) }^ giPijx, t), (3) 

i=l i=l 

where k p is the polymerization rate and D m is the diffusion constant for G-actin. Note that the polymerization 
depends on the total filament density (or number) since each filament within the bundle polymerizes to extend the 
entire bundle. 

To complete the description of the system, the equation for pi (x, t) with i > 1 is 
d t pi{x, t) = -v p Vpi(x, t)+C(x)pi (x, t)+y^ X kti - k pk{x, t)pi- k (x., i)-2 ^ A life pi(x, t)p k (x, t)+(i+l)cp i+ i (x, t)-icpi(x, t). 

k=l k=l 

(4) 

As above, the first term on the right accounts for the ballistic motion, the second term accounts for the angular 
diffusion of the tip. The third term describes those bundles of size and k and i — k coming together to become a 
bundle of size i, while the fourth term describes those bundles that were of size i at the previous time step now 
becoming larger via bundling. We will compute A^fc in the next section within a mean field approach. 

The above equations constitute a general scenario. For the remainder of this work, we will focus on a quasi-two- 
dimensional system — quasi-two-dimensional in that the filaments can overlap one another and so we do not take 
into account excluded volume. In the initial experiments, stars were detected in regions of the solution with more 
spread [l5| . i.e. thinner regions. Hence our reason for studying a quasi-two-dimensional geometry. In this particular 
geometry, we assume that the Arp2/3 initiates a branch in the plane via binding to cither outer edge of a bundle tip, 
hence gt ~ 2. We note, however, that the more recent experiments were conducted in three-dimensional solutions. 
We will address the application of this model to three-dimensions in the discussion section. 
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III. MEAN FIELD THEORY 

Before presenting the results of our kinetic simulations, we introduce a set of simplifications of the two-dimensional 
model to determine the scaling form of the bundle size distribution, ultimately. This result can be tested in the kinetic 
simulations. The first simplification allows us to compute the bundling rates explicitly. A second simplification allows 
for comparison with previous work. A final simplification reduces the model to a single set of coupled differential 
equations. 

The first simplification involves taking the spatially-averaged limit of Eqns. (l)-(4). In mean field, Eqns. (l)-(4) 
become 

^ =2bpl n (t)a(x 7 t)Y / Mt)-cp 1 (t) + 2c P2 (t)-2\ hl pl(t), (5) 

i=l 



^M = -2bpl(t)a(t)J2Pi(t), (6) 
= -k p J2iPi(t)Pm(t) - 2bp 2 m (t)a(t)J2Pi(t), (7) 



dt 



dpi(t) 



, - / j Xk,i-kpk(t)pi-k{t) -2y^\j jk pi(t)pk(t) + + l)cp i+1 (t) - icpijt) for i > 2. (8) 
fc=i fe=i 

A. Tip-tip versus tip-side bundling 

We can now compute A^fc, which is the tip-tip collision rate for two individual ballistically moving filament tips 
in two dimensions. It is given by A*' 1 (t) = (9 c /180°)(2da + dx)v re i(t), where do is a G-actin diameter and dx is a 
crosslinker diameter, v re i{t) is the relative velocity between two growing filaments, and (0 C /18O°) accounts for the small 
angle intersection required for bundling. The relative velocity is ~ p'"[o) v p^ wncrc v p 1S the polymerization velocity. 
Now that we have calculated Af 1; the three-dimensional cross-section for a bundle of size i projected into the plane 
scales as (do + dx)Vi- Therefore, the tip-tip collision/bundling rate for two bundles of size i and size fc coalescing 
is Af fc = ^Ao(t)(\/I+ Vk) with Xo(t) approximated as (6 c /180°)(dG + dx) p™(o) v p- Note that we have neglected the 
cross-linker diameter joining the two bundles, and we have assumed that the cross-linkers can reorient arbitrarily large 
bundles. As for experimentally relevant numbers for X\j, we use d& = 0.003 pm, dx = 0.007/xm, v p = 0.3 pm/sec, 
9 C = 20 degrees (22|, and p m (0) = 10 pM. The dx we have chosen is roughly the size of fascin. The diameter of 
alpha-actinin is approximately 30 nm. Both fascin and alpha-actinin were used in the initial experiments. 

As for our second simplification, we set the concentration of Arp2/3 to be constant. In spatially homogeneous 
system, we cannot study the spatial localization of Arp2/3 presumably required for star formation. Since branches 
also require G-actin to extend, we keep track of the time-dependence of the G-actin so as to observe an upper limit 
on the bundle size distribution, i.e. that the bundle size distribution eventually stops evolving, as is indicated in the 
experiments. We now have 



dpi jt) 
dt 



X li-k{t)pk{t)pi-k{t) -2Y J X?k(t)p l {t) Pk (t) + S hk bpl(t)Y,fr(t) + (i + l)cp*+iW - icpiit) (9) 

k=l k=l i=l 



dp m (t) 
dt 



-kp if* V)Pm (<) " bpl (*) ^Pifat), (10) 



= 1 



where b = 2ba(0) = 2ba. Equations 9 and 10 describe an aggregating system with a time-dependent aggregation 
rate, a time- and concentration-dependent source term, and a size-dependent death rate via capping. As G-actin 
polymerizes into bundles, bundle extension slows. Once the G-actin is used up, the bundles cannot extend any further 
to collide with other growing bundles and the system stops evolving. 
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We study Eqns. (9) and (10) numerically. To do this, we first discuss units. We normalize any density by pi(0) = 
0.002 /uM, the initial density of filaments. This initial density is obtained using a G-actin concentration of 10 \iM and 
a dimerization rate of 3 x 10~ 6 pM~ 1 s~ 1 allowing the solution to dimerize over the course of about 100 seconds. What 
number density does 0.002 pM correspond to? Since we are assuming a two-dimensional geometry for the collision 
process and one pM is roughly equivalent to 600 particles per /im 3 , the effective two-dimensional concentration is 
0.002 pM (600 tips/ pm 3 ) (0.1 pm) = 0.12 tips/pm 2 , where 0.1 /im is the estimated thickness of the star. We use this 
number density when normalizing the two-dimensional collision rate. Now we normalize the time units. Since the 
polymerization of the filaments is the fastest time scale, using k p = 10 pM~ 1 s~ 1 and pi(0) — 0.002 pM, we arrive at 
a time scale of 50 s, which is in rough agreement with the experiments [l6j . 

Using the experimentally relevant parameters listed above, as well as b = 0.1 pAI~ 3 s~ 1 (34|, we numerically integrate 
the mean field system of bundling equations via the Euler method. We impose an upper cutoff, /, on the number of 
bundles and check that the step size of 0.001 is small enough such that pi(t) > as was done in Ref. [Hj]. We have 
checked that our results are essentially independent of the step size and /. With this method, we find that for the 
relevant concentrations of Arp2/3, the tip-tip bundling rate is not large enough to accommodate the excess of single 
filaments, at least for the initial condition of pi(0) = <5;,i in normalized units, the monodisperse initial condition. 
More precisely, the average bundle size decreases from 1.05 to 1.01 as the concentration of Arp2/3 is increased from 
0.005 pM to 0.05 pM. Thus, the tip-tip bundling rate is not large enough as compared to the source rate to induce 
substantive bundling within the fixed timescale of the experiment. 

To arrive at a finite typical bundle size within the time scale of the experiment for the monodisperse initial condition, 
we include another bundling process — tip-side bundling. Such a process occurs and, presumably, it occurs at a higher 
rate then the tip-tip bundling process due to the added length dependence. The tip-side bundling rate for a bundle 
of size i coalescing with another bundle of size k is given 

by K\ = Ao(t)(V* + where L(t) is average length 
of the "side" bundle. In units of the monomer size, the equation for L(t) is 

dL 

-ijT = k p p m {t). (11) 

As for -L(O), the average length of the initial seeds is not necessarily well-characterized. The dimerization rate used 
above assumes that L(0) = 2, therefore, that is what we implement. 

For tip-tip bundling, the average bundle size decreases with increasing <r(0) due to the uptake of G-actin into 
branching. However, for tip-side bundling there is an added effect. Namely, since the length of filaments is dependent 
on the concentration of G-actin, the more Arp2/3 there is, the smaller the average length of filaments, the smaller the 
tip-side bundling rate resulting in a smaller bundle size typically. This dependence is in accordance with a molecular 
dynamics simulation based on the later star experiments [17j . Certainly, this trend is observed for large Arp2/3 
concentration where stars no longer form since the filaments remain short and only a dendritic cloud forms. However, 
in the experiments, it has been observed that a critical Arp2/3 concentration is needed to initiate star formation. 
This is presumably because there has to be enough filaments around to bundle. In our mean field model, this latter 
dependence can be accounted for by decreasing the initial filament concentration to such a level that the typical time 
scale becomes slower and slower such that depolymcrization has to be taken into account. Since we do not take 
depolymerization into effect, we presumably we will not be able to account for a threshold concentration of Arp2/3 
in bundle formation. 

The length dependence of the tip-side bundling rate is not as straight forward as above. For small angle collisions 
to occur, the entire length of a filament/bundle is not necessarily available for bundling. To estimate the effective 
length for bundling requires knowledge of the spatial distribution of the initial individual filaments and of the Arp2/3. 
In the absence of Arp2/3, for a distribution of initial single filaments whose radial position is much smaller than the 
final typical filament length, tip-tip bundling is more likely to occur than tip-side bundling. However, for a more 
broadened distribution of initial seeds, the ratio of tip-side to tip-tip collisions will increase. The spatial distribution 
of the Arp2/3 must also be taken into account since the more broad the distribution of Arp2/3, the ratio of tip-side 
to tip-tip collisions may potentially increase. Kraikivski and collaborators [24[ incorporate this effect into the tip-side 
bundling rate by decreasing the critical bundling angle. Since spatial information is crucial for estimating the available 
length for bundling, we merely study pi as a function of 8 C . 

Now that we have fully defined the various parameters in mean field, we integrate Eqs. 9 and 10 with tip-side 
bundlingand neglect tip-tip with c = for the remainder of the section. When we take into account tip-side bundling, 
we obtain average bundle sizes larger than unity within the time scale of the system. 

In Figure [21 we plot the final bundle size distribution for different concentrations of Arp2/3. To plot the bundle 
size distribution, or the probability of having a bundle of a particular size, we plot tPz ip . as a function of i. The 
denominator is needed for normalization since the mass density is not conserved. Plotting the final bundle size 
distribution on a log- linear plot, we see that eventually the distribution falls off exponentially with bundle size with 
the decay depending on the concentration of Arp2/3 (sec Fig. It is this decay that sets the average length. 
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We also plot the average bundle size as a function of Arp2/3 concentration (Fig. 3(a) |. As expected in mean field 
with no depolymerization, an increase in Arp2/3 concentration decreases the average bundle size due to the faster 
decay of G-actin. In addition, the faster G-actin gets soaked up into branches, the smaller the average length of 
filaments/bundles, thereby, again, decreasing the bundling rate and, hence, the average bundle size. We also note 
that bundles are able to form in the absence of Arp2/3 with this initial pi(0). 

In Figure [3(b) | we plot the average bundle size as a function of 8 C . This parameter is an attribute of the cross- linker 
and so can be studied by modifying/changing the cross- linker. The bundling rate is proportional to 6 C . so an increase 
in 6 C increases the average bundle size. However, the dependence is complicated by the fact that the source term 
depends on the total number of bundles such that the more bundling, the less branching, the longer the filaments 
thereby increasing the average bundle size even more so. 

Comparison with previous work: To make a more direct comparison with the work of Kraikivski and collabora- 
tors [Uli we simplify the branching/source term to bp m (t) Pi (0) = bp m (t), where b = 2bap m (0). With this 
simplification, we obtain exact expressions for p m (t), L(t), and M(t), the total mass density. In particular, we have 



Pm(t) = ^y^(scch(-fc pV /2A(i + D 2 )))-' 2 (12) 



M(t) = y/2D' 1 tanh(h p y/2D' 1 (t + D 2 ))-^, (13) 

kp 

with Di and D 2 determined by initial conditions and k p — k p p\{Q). Once we know M(t), we also know L(t). Replacing 
L(t) with its long-time value, we analyze 

= A (i)L(oo) J2(Vk + x/T - k)p k (t) Pi . k (t) - 2A (i)L(oo) ^(Vfc + V~i) Pl (t) Pk (t) + S iA bp m (t). (14) 
fe=i fc=i 



| = Ao(t), d5) 



Rcparameterizing time with 
we arrive at 

= L(oo)(J2(Vk + sfi=k)p k { T ) Pi - k {T) - 2^(v / fc + Vi)pi(T)p k (r)) + 5 lA b> (16) 

k=l fe=l 

with b' oc b. 

Numerical simulations of this further simplified model show that tip-tip bundling alone is not enough to produce 
substantive bundle sizes and so we do not expect the tip-tip bundling to be the dominant mechanism of bundling. 
Kraikivski and collaborators [24| reach a different conclusion. We suspect this difference is due to the fact that our time 
scale involves the polymerization rate, the tip-branching rate, and the initial concentration of individual filaments. 
The time scale computed by Kraikivski and collaborators [24J does not appear to involve the initial concentration 
of individual filaments. However, since Arp2/3 must branch on pre-existing filaments, we expect the time scale to 
depend on this parameter. To arrive at a time scale comparable to the experiments, the branching rate for Kraikivski 
and collaborators is several orders of magnitude smaller than ours such that their branching/ source rate is more 
comparable to their tip-tip rate so larger bundles are able to form in their model. 



B. Scaling form for bundle size distribution 



The numerical results of the previous subsection suggest an exponential form for the bundle size distribution, at 
least for higher concentrations of Arp2/3. Interestingly, an exponential tail has been measured for an actin-fascin 
system [26|. While a finite of amount of G-actin accounts provides for an upper bound on the time evolution of 
the system and, therefore, determines the radial extent of the system, we expect the bundle size distribution to stop 
evolving before all of the G-actin polymerize (in the absence of zippcring). Therefore, we approximate the G-actin 
concentration to be constant. Such a system can be constructed by supplying the homogeneous system with G-actin 
at a rate that equals the loss of G-actin to polymerization and branching. 
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If we assume a constant G-actin concentration, then L(t) increases linearly in time. We are free to reparamctcrize 
time using 

§ = Xo(t)L(t), (17) 
with Ao (t) normalized appropriately such that r ~ t 2 . Then equation for the density of bundles sizes becomes 

= J2(Vk + V^ r k) Pk (r)p i - k (T) - 2]T(7fc + Vt)p t (r)p k (r) + 4,i 5 (r), (18) 
fe=i fc=i 

with g(r) = ^ Prn ^ N(t), where N(r) = J2iPi( T )- A simpler form of this equation was studied by Davies and 

y p Ap ( 0)t 

collaborator J3a] where the time-dependent source term did not depend on the density. We extend those results here. 

Before doing so, in the absence of the source term (ns), there is no exact solution to date. However, since the 
bundle size dependence of the collision rate is homogeneous of deg ree p, i.e. V~i + Vk = k p (1 + (k/i) 1 Z 2 ) with p = 1/2, 
there has been much study of the scaling form for pi (r) [U |37| The dynamic scaling hypothesis [2l|, [37j predicts 
that 

Pi 3 (r) ~ (19) 
s(t) z s(t) 

at large times with s(t) ~ r 2 . The s~ 2 prefactor enforces mass conservation and the scaling of s with r can be derived 
by using Eq. 19 in the continuum limit of the kinetic aggregation equation given by 

ry roc 
d T p{y,r) = dz(z 1 / 2 + (y-z) 1 ^)p(z,r)p(y-z 1 T)-2 dz{z l l 2 + y 1 ' 2 )p{y,T)p{z,r). (20) 
Jo Jo 

Moreover, ^S(x) is expected to fall off at least exponentially at large x, and for small x, ^(x) ~ .t~ 7 with 1 < 7 < 
3/2 [38|. The lower bound enforces the normalizability of the distribution. The upper bound comes from asymptotic 
analysis of the kinetic aggregation equation using the dynamic scaling hypothesis in the limit of small x. Tighter 
bounds can be computed [3g|. We refer the reader to the excellent review by Leyvraz [2l[ where this analysis is 
presented in detail. 

And now for the analysis in the presence of branching. For small bundle sizes and long enough times, one would 
expect the aggregation terms and source term to balance one another. In other words, one would expect a pscudo 
steady state (pss) to emerge with 

pT S - \W)r A , (21) 

with A = 7/4, which is obtained from the constant source case 39]. For intermediate bundle sizes and long enough 
times, one expects a self-similar solution of the form 

Pi (r) ~ r- 7 /^( i / T i/2 )) (22) 

where we used (1) the conservation of mass in the presence of a time-dependent source, or Yli=i iPi( T ) ~ 
J drN^/r 1 / 2 , and (2) the scaling form for Pi(r) consistent with the source- free aggregation equation with this 
particular form of the collision rate, or Pi(r) ~ r~( 1+ ^ 9 ^(i/r 9 ). For this intermediate range of i, a scaling collapse is 
demonstrated in Figure [4 ( a) | With this scaling for pi(r) in the intermediate region, the integrated newly introduced 
mass eventually deca ys wit h time and the system approaches the constant mass (no source) case with Pi(r) scaling 
as /9™ s (r). See Figure [4(b)] for a scaling collapse for this large bundle size region solution. 

While we may not observe all of these regimes in finite-dimensional simulations of such a system, since presumably 
there is an upper bound on the bundle size, based on the above analysis, it is reasonable to expect that 

-A 



Pi ~ i exp(—i/ich), (23) 



where i c h is some characteristic bundle size that depends on the amount of Arp2/3, for example, and A « 7/4. 
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IV. KINETIC SIMULATIONS 

Now that we have a prediction for the scaling form for the bundle size distribution in the spatially homogeneous 
system, let us test what aspects of that prediction will persist in the kinetic simulations with spatial degrees of 
freedom. We will also measure the local orientation of the filaments to determine if indeed a star has formed as 
opposed to a network of bundles, which arc not diffcrcniated in the analysis above. To accomplish these tasks, we 
perform kinetic aggregation simulations in quasi-two-dimensions. More specifically, we allow for a third dimension 
only in that we allow growing filaments/bundles to cross each other. Furthermore, we not only keep track of small 
angle cross-links that leading to bundling but also large angle cross-links as well. These large angle cross-linkers do 
not reorient the filaments, but account for connected clusters of filaments. These connected clusters will become 
important when including the finite-size of a bead embedded in the system to closely model the initial experiments. 
We do not incorporate dynamics of the linkers. Once two filaments are bound, they remain bound for the time scale 
of the experiment. 

A. Algorithm 

To simulate the actin-Arp2/3-cross-linker system we keep explicitly keep track of the filaments and the Arp2/3, but 
not the G-actin and the cross-linkers. In other words, the concentration of G-actin and cross-linkers are parameters, 
while the positions of the actin in filaments are updated and recorded as are the positions of the Arp2/3. Explicit 
tracking of the Arp2/3 allows us to test the assumption that localization of the Arp2/3 is required for star formation 
since Arp2/3 promotes dendritic growth as opposed to oriented, bundled growth. As for the constant concentration 
of G-actin, in a spatially homogeneous system, the filaments bundle indefinitely. However, in a system with spatial 
degrees of freedom, the bundling rate acquires spatial dependence such that as bundles orient as they extend outward 
from the center of the system, the bundling decreases towards the perphiery of the system in even the presence 
of G-actin. We expect this spatial dependence, and not the G-actin concentration, to limit the bundle size in 
the absence of zippering. Of course, a finite amount of G-actin is required to limit the radial extent of the star. 
Since previous numerical studies have focused on this aspect, we will not focus on it here fl7j and assume that 
filaments grow at the same rate near the center of the system. As for the constant concentration of cross-linkers, 
the initial experiments found cross-linkers both the dendritic and bundled regions thereby justifying our simplification. 

We perform the following steps: 

(1) Initialize Np single filaments at some random position with the plus end pointing at some random angle uniformly 
distributed between [0, 2tt]. The random positions are chosen uniformly within an annulus whose inner radius is i?t>, 
the radius of the bead, and whose outer radius is Rb + Ri- 

(2) At each simulation time step, a filament grows by 100 nanometers, or approximately 30 monomers. This 
compression of data storage will allow us to keep track of percolation (connected cluster) properties of the cross- 
linked network over length scales of the experiment. In addition, at each simulation time step, the angle of the 
segment is varied by a Gaussian random angle with a variance of 3 degrees with respect to the angle at the previous 
simulation time step to approximate a persistence length of tens of microns. When i?f, > 0, a filament cannot 
grow in the bead region unless the entire filament and the filaments connected to it via small-angle and large-angle 
cross-links, the entire cluster, can be pushed backward/outward away from the bead. In other words, the entire 
connected cluster ratchets off the bead. To keep track of the clusters of cross-linked filaments, small angle and large 
angle, we implement the Hoshen-Kopelman algorithm ■ 

(3) After updating the growth of each filament, a check is done to determine if any two new segments of two 
growing filaments intersect at an angle less than 9 C or are within a distance dx + da of another filament and 
are about to intersect at an angle of less than C . If so, then the filament in the smaller bundle is reoriented 
towards the filament in the largest bundle. A check is also made if the new segment intersects with any older 
segment of another filament. If so, the new segment realigns with the older segment (provided the intersection 
angle is less than 9 C ) and a bundle is formed. For both checks, new with new and new with old, if the intersection 
angle is greater than the two filaments become "permanently" connected via a cross-linker, but no realignment is made. 

(4) As for the Arp2/3, the position of Na Arp2/3 is initially randomized in the same way that the initial fil- 
ament seeds are. The Arp2/3 subsequently diffuses with a diffusion constant corresponding to approximately 1 
fim 2 1 s. After updating the positions of the Arp2/3 in parallel, a check is made to determine whether an Arp2/3 is 
within a distance dA of a filament. If so, a new filament is nucleated 70 degrees from the growing end of the filament 
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and the position of the Arp2/3 is recorded. 



(5) Steps (2) through (4) are repeated until all the Arp2/3 binds. Then, the positions of the growing ends of 
the filaments are recorded. Also, the local orientation of the filaments as a function of distance from the center is 
measured via the following method: 

Two segments are chosen at random, if the two segments are within a distance df of each other, then the cosine of 
the angle between the two segments is calculated and binned according to the distance from the center. 

(6) Since we are interested in how bundles at the extremities of the system emerge from the center of the 
system, those filaments whose growing end point outwards and exceed a radius R cu t and are tagged and continue to 
grow via (l)-(4) until the bundle size distribution, which is computed at each time step, stops evolving. The local 
orientation of the filaments is measured again further out. A final bundle size distribution is also outputted. The 
final bundle size distribution will depend on R cu t- We will use this dependence to say something about the nucleation 
of bundles from the center of the system. 



(7) (l)-(6) is repeated 1000 times for sample averaging. 



Two snapshots of the simulations, one without a bead and one with a bead, are depicted in Figs. 5(a) and 

m\ 

We also note for computational simplicity we do not explicitly increase the bundle size when filaments bundle. The 
reorientation of the individual filaments for bundle lengths smaller than the persistence length is sufficient to keep 
the two filaments aligned. In other words, we did not capture the y/i bundling rate dependence explicitly, however, 
the bigger the bundle, the more likely it is to collide with another bundle due to the increase in local filament density. 
For (6), we also decrease the amount of angular diffusion of the tips by 1/i, where i is the bundle size. This is to 
capture the appropriate persistence length dependence for bundles of filaments. 



TABLE I: Table of Parameters 



Parameters 


Variables 


Values 


Polymerization rate 




10 pM^s' 1 


Branching rate 


b 


0.1 (jM- 3 s _1 


Capping rate 


c 


0- Is" 1 


G-actin concentration (when fixed) 


Pm 


~ 10 pM 


Arp2/3 concentration (when fixed) 


a 


0-10 nM 


Diffusion constant for Arp2/3 




1 pm 2 s~ 1 


Critical bundling angle 


0c 


0-30° 


Diameter of G-actin 


da 


3 nm 


Diameter of cross-linker 


tlx 


7 — 17 nm 


Arp2/3- filament binding distance 


(I a 


0.01 pm 


Local orientation distance 


df 


0.01 - 0.03 pm 


Radius of the bead 


Rt 


0.5 pm 


Max. rad. of initial filament seeds from center - 


Rb Ri 


0.4 — 1.5 pm 


Radius for bundle size determination 


Rcut 


0.75 - 0.95 pm 


Number of initial filament seeds 


N F 


100 - 500 


Number of Arp2/3 


Na 


250 - 650 



B. Results 



1. Absence of bead: Rb = 

Bundle size distribution : 

In Fig. [5] we plot the logarithm of the bundle size distribution for 8 C = 30°. As for the other parameters, Ri — 
0.6 pm, Np = 500 and Na = 250. And, unless otherwise specified, dx — 0.017 pm, do — 0.003 pm, dA = 0.01 pm, 
df = 0.03 pm, and R cut — Rb + 0.75 pm. Note that dx is slightly larger than in mean field, at least for the bulk of 



10 



the simulations. Also note that the concentration of initial seeds is several orders of magnitude higher than in mean 
field depending on the thickness of the star. 

From Fig.[6]find that the bundle size distribution approaches a limit as time proceeds. In other words, the collision 
rate decays substantially with time even in the presence of G-actin over the lengthscale of approximately 5 microns. 
This is because once the bundles form near the center of the system, collect together and polymerize outward in an 
oriented fashion, the spacing between bundles increases, therefore, the collision rate decreases. Bundles are less likely 
to merge at the system's edge. 

We investigate the effect of 9 C on bundle formation. If all other parameters are known, measurement of the bundle 
size distribution allows one to determine the amount of reorientation by the cross-linker — a microscopic detail. From 
the mean field kinetic aggregation calculations in Section IIIB, the bundle size distribution is exponential with an 
algebraic correction. See Eq. 23. In quasi-two dimensions, at least for 6 C > 10°, the data suggests a similar form, 
but the algebraic correction is a power law with exponent approximately A — 1 w 0.6 as opposed to A — 1 = 0.75 
computed in Sec. IIIB. (The bundle size distribution contains an extra factor of the bundle size i in comparison to 
Pi). To see this algebraic correction more clearly, we also plot the log-log of Fig. 7(a) in Fig. |7(b)| The exponent 
obtained in two-dimensions may be in slightly better agreement with the one found for a constant aggregation kernel 
with a constant source term [21|, which is A — 1 = 0.5. This is reasonable given that we have not explicitly included 
the size of the bundle in the simulation. Work by Connaughton and collaborators [4l[ have investigated the role of 
spatial fluctuations for a system of diffusing, coalescing particles in the presence of a constant source term and found 
that the exponent depends on the spatial dimension for d <= 2. It would be interesting to extend this work to the 
semiflexible case to more accurately determine analytically the algebraic correction in finite-dimensions. Furthermore, 
it seems that the spatial distribution of the Arp2/3, which was obviously not accounted for in the mean field theory, 
does not alter the form of the bundle size distribution from mean field. 

We have also measured the bundle size distribution as a function of concentration of Arp2/3 by changing Na but 
keeping Np + Na fixed. Svitkina and collaborators [HI argue that Arp2/3 provides enough filaments for a broad, 
dendritic base from which filaments can emerge and bundle to form larger bundles than in the absence of Arp2/3. 
This concept of a larger base leading to larger bundles is a natural one and is supported by the fact that while it is 
less likely to have smaller bundles the larger the Arp2/3 concentration, the more likely it is to have larger bundles 
thus increasing the average bundle size overall. See Fig. [8] The average bundle size is 15.50 ± 0.32 for Na = 450 and 
18.46 ± 0.40 for Na = 650. In contrast with mean field, here, an increase in the concentration of Arp2/3 increases 
the average bundle size. 

We observe a larger base leading to larger bundles concept when we study the bundle size distribution as a function 
R cu t- The larger R cu t, the smaller the base from which filaments can emerge and bundle resulting in smaller bundles. 
For example, for R cut = 0.85 pm (as compared with R cut = 0.75 fxm), the average bundle size is 9.03 ± 0.18 and for 
Rcut = 0.95 fj,m, the average bundle size is 5.57 ± 0.13. 

While we find support for the convergent elongation model here, at least in the absence of depolymerization, we do 
not find a threshold concentration of Arp2 / 3 for bundling. In addition, we find that in order to achieve bundles of 
reasonable size (of order 10 on average), one needs enough filaments, branched or otherwise, within a small enough 
region. If we increase Ri to 1.8 microns from 0.6 microns, then the average bundle size is 7.58 ± 0.14 as compared to 
14.94±0.20 (with some nontrivial density dependence in between). Beyond Ri = 1.8 microns, the average bundle size 
continues to decrease. Indeed, there is a threshold filament concentration needed for substantive bundle sizes. This 
threshold concentration depends on somewhat on the ratio of Na/Nf but more strongly on Ri. If depolymerization 
is included, then we should observe a threshold concentration of Arp2/3 for substantive bundling. 



Orientation of filaments: We also measure the local orientation of the filaments with respect to one another to 
determine if indeed a star has formed as opposed to an isotropic network of bundles. First, let us address several 
scenarios. For an isotropic system of nonpolymcrizing filaments, in the absence of bundling, the average of the cosine 
of the angle between nearby filaments, 0, should vanish. For a static system of isotropically distributed bundles, for 
length scales on the order of the bundle cross-section, the average of the cosine of the angle between nearby filaments 
should be unity, but on larger length scales of order the bundle cross-section, will decrease. For a static system of 
aligned bundles, will remain close to unity even as the length scale for nearby filaments is increased. For a static 
system of 70° branched filaments, then = 0.3, provided the mesh size of the system is not too small. Haviv and 
collaborators (l6| observe this trend in simulations modeling their experiments. 

What about systems with polyermization? For a system of polymerizing filaments, the filaments that grow into 
new regions of the system tend to be oriented in the same way, even in the absence of bundling. Otherwise, they 
would not be able to reach the same new regions in the first place. Therefore, > for the outer exterminities of 
the system irregardless of bundling. 

So, for star formation to occur, we expect 0(r) to inrease from zero near the center of the system. Given the 
cylindrically symmetric initial conditions, it is more probable for filaments of opposite orientation to overlap near the 
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center. As r increases, O should increase to approximately 0.3 in the dendritic region with this number being mesh 
size dependent. Beyond the dendritic region, towards the extermities of the system, O will contine to increase beyond 
the value of for the 9 C = 0° case. 

Before presenting for the same parameters as in Fig. [JJ we plot 0(r) for a system where we do not expect to 
observe star formation in Fig. [9] More precisely, Np = 100, Na = 650, Ri = 2 /jm and 9 C = 0. We expect this system 
to Arp2/3 dominated such that is near 0.3. The spatial distribution of the Arp2/3, denoted by Prob a , is also 
plotted. We see that the average local filament orientation modulates around 0.3. If the Arp2/3 is more localized, 
then we expect the average local orientation between filaments beyond the Arp2/3 region to increase beyond 0.3. To 
check for this, we let Ri = 1.5 fim and compare with the Ri = 2.0 /im results. We sec that beyond 1.5 microns, the 
Arp2/3 decays and the average local orientation between filaments increases beyond the Ri = 2.0 /im data. 

Now we plot 0(r) as a function of 9 C for the data depicted in Fig.[7J See Fig. [TU] For 9 C = 0°, where no bundles are 
formed, 0(r) increases from zero as r is increased. We do not observe a plateau at 0(r) = 0.3 in the Arp2/3 region 
(shown later). Even if we increase Na relative to Np, due to the high concentration of filaments with df = 0.03 fim, 
we do not observe the plateau. Beyond the Arp2/3 region, continues to increase, until ultimately, beyond 2 microns, 
begins to decrease indicating that filaments begin to overlap again as they grow. 

How does depend on 6 C > 0°? For df = 0.03 [im, it turns out that is slightly smaller for 6 C = 30° in comparison 
to C = 0° near the center of the system. This result is a bit counterintuitive, but reasonable. The reorientation 
of filaments biases the system towards a network of bundles as opposed to a branch dominated system in this high 
density regime. Since df is larger than the mesh size, we are more likely to measure filaments/bundles of different 
orientation for 9 C = 30° than for 9 C = 0°. Decreasing df increases for 8 C = 30° near the center of the system such 
that the 9 C = 0° and 6 C = 30° curves approach each other. Near the periphery of the system, as 8 C increases from 
zero, the asymptotic value of increases. The asymptote of approximately 0.95 for 9 C = 30° certainly indicates star 
formation when compared to the C = 0° data. 

Summary: In the absence of a bead, 9 C must be sufficiently large enough to increase the bundling rate to 
approach an average bundle size of approximately 10 filaments per bundle. Furthermore, there is a threshold 
concentration of filaments, branched or unbranched, required for substantive bundling. There exists a spatial 
contribution to the collision rate for bundling that does not exist in mean field (by construction). In addition, 
for star formation to occur, as opposed to an isotropic network of bundles, Arp2/3 needs to be localized with 
bundles/filaments polymerizing beyond the Arp2/3 region. Only in regions where the Arp2/3 is not present can 
can attain a value near unity. Also, reorientation is needed for to approach unity. 

2. Presence of bead, Ri = 0.5 /im 

Now we investigate the effect of a finite-sized bead on the bundle size distribution and the average local filament 
orientation. The bead not only acts as a steric object excluding filaments, but growing filaments can "ratchet off" the 
bead. More specifically, a filament with its growing end toward the bead fluctuates away from the bead such that a 
monomer can be added. Once the monomer is added and the filament fluctuates "back" toward the bead, the entire 
filament moves backward or away from the bead. If this filament is connected to other filaments, then those filaments 
move backward by the same amount as well. So, this ratcheting effect will expand the dendritic region, at least under 
certain conditions. 

Bundle size distribution: First, we investigate the bundle size distribution. As in the Rb = case, the bundle 
size distribution stops evolving substantially after the bundles have localized enough to decrease the collision rate 
substantially towards the extremities of the system. See Fig. [IT] 

We also study the bundle size distribution as a function of 9 C in the presence of the obstacle. See Fig. [12j As in the 
no bead case, the larger 9 C , the larger the average bundle size. More specifically, for 9 C = 10° with Ri = 1.2 ^m, the 
average bundle size is 1.55±0.01, for 9 C = 20°, it is 3.48±0.01, and for 9 C = 30°, we measure 7.20±0.08. If we choose 
Ri = OAfim, then the initial densities of filament seeds are approximately the same for Rb = and Ri = 0.6 fim 
case (3.6-7T /im 2 versus 3.9-7T /xm 2 ). Then, we have for 9 C = 30°, an average bundle size of 11.23 ± 0.27. The average 
bundle size is smaller when compared to the system with no bead as expected. The obstacle of the bead makes it 
more difficult for filament on the other side of the system to grow and eventually meet up with growing filaments of 
the same orientation. 

The bundle size distribution data suggests that the form of the bundle size distribution is, as in the Rb = case, 
a decaying exponential with a power law correction. For 9 C > 10°, the exponent A — 1 = 1.0 in comparison to the 
Rb = case for same initial density of seeds. This power law correction depends on the density of initial filament 
seeds. By increasing Ri from 0.4 microns to 1.2 microns we find that the exponent of the power-law correction varies 
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from approximately A — 1 = 0.8 to A — 1 = 1.4 for Ri — 1.2 \im. 

It turns out that the average bundle size decreases with increasing Arp2/3 concentration at least for a range of 
Ri near 1.2 fim. This result differs from the Rb = case. However, the observation of a threshold concentration of 
filaments still holds. The bundle size distribution also depends on R cu t as in the Rb = case. The larger R cu t the 
smaller the base of bundle and the smaller the average bundle size. Finally, we also decrease dx to 0.007 [im and 
find that the average bundle size decreases to 6.76 ± 0.06. 



Orientation of filaments: Now we investigate 0(r) in the presence of the bead. See Fig. 13(a) Since the filaments can 



no longer access the center of the system, 0(r = Rb) > 0. For r > R b , the ratcheting effect extends the dendritic phase 
of the filaments such that 0(r) is nonmonotonic for r < 2 microns. This nonmonotonicity is due to the ratcheting 
effect where branched filaments get pushed outward such that there is some "memory" of the randomized initial 
conditions further away from the bead then in the Rb = case. See Fig. [10] for comparison. 

For 0(r) even further away from the bead, r > 2/im, we notice a second nonmonotonicity (as compared to the 
Rb = case). If we plot as a function of C , we notice that the negative slope for the edge filaments increases in 
magnitude with 9 C . Due to the ratcheting mechanism, filaments are less ordered emerging from the isotropic region 
such that there is a larger spread in orientation among the bundles (and individual filaments). This large spread in 
orientation from bundle to bundle eventually leads to overlap between bundles at an angle not necessarily less than 
9 C . Hence, one observes a decrease given the overlap between different bundles. Moreover, the obstacle of the bead 
prevents some growing filaments from "meeting up with each other" even if their seeds are close to each other and 
have a similar orientation. So 0(r) is smaller for Rb > as compared to Rb = 0. If we increase Ri, we decrease the 
ratcheting effect and the become less nonmonotonic. See Fig. |13(b)| 

Capping: Finally, the initial in vitro experiments investigated the effects of capping protein on the star forma- 
tion. Capping protein needed to be depleted by a sufficient amount in order to observe the star formation. However, 
an increase in the concentration of capping protein, for small concentrations, induced further cloud growth, i.e. 
the radial extent of the dendritic/cloud increased [151 ] . This result is counterintuitive because capping inhibits 
polymerization so one would think that the radial extent of the network would decrease. Eventually, for high enough 
capping protein concentration, the radial extent of the network decreases such that even the formation of clouds is 
inhibited. So, the enhancement of cloud growth is not so obvious. How do we explain this? 

The interplay between the ratcheting effect of the bead and the percolation of the cross-linked filaments/bundles can 
explain this result. For small capping, filaments are longer and, therefore, have a great chance of forming large-angle 
cross-links with other filaments. Increasing large-angle cross- links increases the likelihood that filaments/bundles 
are connected to other filaments on the opposing side of the bead, i.e. percolation across the system is more likely. 
Then, shorter filaments growing into the bead cannot ratchet backward because the filaments on the opposing side 
of the bead cannot move into the bead. However, as capping is increased, there is less connectivity and, hence, less 
chance of connectivity between two sides of the bead, thereby among allowing the ratcheting effect to dominate and 
enhancing the radial extent of the network. We observe an extreme case of this in Fig. [TU Two clouds have been 
"born" if you will, for a capping rate of 0.6 s _1 . See Fig. [I])] for the averaged data. As capping is increased from 
zero, the position of the filament's growing tips becomes more extended. Increasing Arp2/3 will mimic this effect, 
provided the filament density is not too large, since filaments are shorter on average thereby decreasing the likelihood 
of forming large-angle cross-links. 

Summary: For Rb > 0, again, we find that the bundling rate and the filament density needs to be sufficient 



enough to attain average bundle sizes of order 10. The bundle size distribution has the same form qualitatively as the 
Rb = case. While the bead does inhibit bundle formation to some extent, less ordered star formation still occurs 
with smaller than the Rb = case. The presence of the bead also leads to a nonmonotonic capping protein effect 
via the interplay of percolation and ratcheting. 



3. Ordered initial conditions and Rb = Na = 



Finally, we go beyond modeling the experiments and we investigate the case when bending of the actin filaments 
over long enough lengthscales determines the bundling rate. This case can be studied with seeds emerging from an 
ordered array with the same orientation. At long enough lengthscales, the filaments will ultimately bend and collide 
with one another at small angles and form bundles. See Appendix A for analysis of the bundle size distribution in 
this case. The simplified analysis yields a bundle size distribution that is a nonmonotonic in the bundle size, which 
differs from the disordered initial condition and small length scales case. However, if we incorporate the effect of 
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reorienting filaments via bending into the disordered initial condition model, we also find a nonmonotonic bundle size 
distribution with bundle size. 

We test this nonmonotonic prediction numerically for Np — 100. See Fig. [TH The seeds start 0.1 microns apart 
and dx = 0.017 lira. As 8 C increases to 30°, for filament lengths of 8 mm, the bundle size distribution starts to become 
nonmonotonic. In other words, the distribution begins to look qualitatively different from Na > 0, randomized initial 
conditions and shorter length scales. We expect that such a distribution can be obtained experimentally as long as 
the concentration of G-actin can maintained over scales much larger than the persistence length. 

V. DISCUSSION 

We have extended the kinetic aggregation approach to study the morphology of an actin-bundling-protein-Arp2/3 
system. Such an approach began with Smoluchowski [2(| and has been used to describe numerous systems ranging 
from coagulating colloids to polymers to ordinary and galactic dust to name a few (2lj . Kinetic aggregation modeling 
is orthogonal to the equilibrium approaches of actin bundling via electrostatics (42j . Mechanisms for setting a typical 
bundle size in these scenarios range from the long-range electrostatic repulsion between like-charged filaments [44| 
to the interplay between electrostatics and steric effects [43| to a energetic penalty for chiral packing [45|. However, 
these ideas are applied to solutions of fixed length actin filaments with no branching. The star system contains both 
polymerization and branching, suggesting a nonequilibrium approach. 

With minimal assumptions (aside from assuming the system is kinetically driven), we have computed the bundle 
size distribution for a system of polymerizing actin filaments that branch and bundle. In mean field, in contrast with 
Kraikivski and collaborators [24| , we find that tip-tip bundling alone is not enough to generate bundles of substantive 
size. Moreover, in our quasi-two-dimensional simulations, we keep track of tip-tip (new-new) bundling versus tip-side 
(new-old) bundling and find for most initial conditions studied that the tip-tip bundling events are approximately 
four to fives lower than the tip-side events. While not as dramatic a difference as in mean field, there is indeed a 
difference. 

In finite-dimensions, the average bundle size is controlled by a spatial factor in the collision rate, which limits 
the bundle sizes in a nontrivial way as the system extends itself. Since the limiting factor comes into play near the 
extermities of the system, the notion of the lambda percursor is an important one for substantive bundling to occur. 
Bundles should form early on within the dendritic meshwork and extend themselves as well as collide with nearby 
bundles soon after emerging from the dendritic network. In other words, our results support the convergent elongation 
model. 

The conjecture that Arp2/3 is crucial for filopodia formation is not bourne out in these results, while a threshold 
concentration of filaments, branched or otherwise, is, at least in the absence of depolymerization. In recent studies of 
neuronal cells, it has been shown that the depletion of Arp2/3 results not only in less lamellipodia but in less filopodia 
as well [5l| . again, highlighting the importance of Arp2/3. Extending our simulation to include depolymerization as 
well as excluded volume interactions would be useful for further studies of the conjectured importance of Arp2/3. 
The excluded volume interactions may help explain how Arp2/3 localizes, even in the absence of the bead, to form 
stars. One could also take a step back from the experiments and explore the role the branching angle of the Arp2/3 
on filopodia formation. Edelstein and Ermentrout (4r| have investigated the interplay between the branching angle 
and cross-linking and found that for small branching angles, an instability in the isotropic state emerges such that 
filaments orient along one direction and so are not able to bundle effectively at least for short filaments. One may 
expect the larger the branching angle, the larger the base of the lambda precursor, resulting in larger bundles. 

As for the shape of the bundle size distribution, we determine in mean field that the distribution is exponentially 
decaying but with a power-law correction. The crossover between the two regimes can be tuned by varying the 
concentration of Arp2/3. The quasi- two-dimensional simulations exhibit this form as well but the exponent depends 
on various factors such as the presence of the bead. Moreover, in the quasi-two-dimensional simulations, increasing 
the concentration of Arp2/3 (keeping the G-actin concentration fixed) increases the potential base for the lambda 
precursors, at least in the absence of the bead, leading to larger bundle sizes on average. We begin to obtain a 
qualitatively different bundle size distribution for filaments emerging from an ordered array, where filaments can only 
collide via transverse bending. The bundle size distribution becomes nonmonotonic in the bundle size for a large 
enough bundling rate. Nonmonotinicity can also be accounted for by introducing a cutoff in the number of filaments 
the cross-linkers can reorient given the cost in bending energy interestingly enough. Determining whether the bundle 
size distribution is peaked or not should presumably be the next experimental test of the star system. We should 
note bundle size distribution measurements of a network of bundles in an actin-fascin system suggest that there is a 
peak in the distribution of bundle sizes, though evidence for the peak is essentially one data point [26| . 

Initial experiments demonstrated that star formation was more likely to occur in places where there was more 
spreading of the solution across the coverslip, hence the quasi-two-dimensional modeling. However, subsequent ex- 
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pcriments and simulations have studied star formation in three-dimensions [la . Il7| . The model studied here can be 
extended to three-dimensions. In the mean field theory, the parameter introduced in Equation 1 will now be size 
dependent since the Arp2/3 could branch off of the surface area of the bundle as opposed to the two outer filaments 
in each bundle in two-dimensions. The larger the bundle, the more likely for Arp2/3 to branch off its surface area. 
This increase in the nucleation of new filaments will offset, to some extent, the decrease in the collision rates in three 
dimensions and presumably account for star formation in three-dimensions as is observed. In three dimensions, one 
expects tip-side bundling to dominate over tip-tip bundling to an even greater extent. 

While we have focused on bundle formation via kinetic aggregation in the absence of a membrane, others have 
begun to study the interplay between bundles and the membrane. Atiligan and collaborators [47| propose a filament- 
membrane interaction and show that filopodia extension rate is faster for an elastic membrane as compared to a rigid 
membrane. Lan and Papoian |48| build on the previous simulations by incorporating monomeric flow to demonstrate 
that rate of filopodia extension is very sensitive to the retrograde flow of monomers. Other work has investigated the 
possibility of the membrane inducing bundling j50| . It would be interesting to extend our work to include a dynamic 
obstacle and determine its effect on bundle formation. We should note that a more recent simulation by Zhuravlev 
and Papoian studies in detail the effects of the binding and unbinding of capping protein on filopodia formation 
and discovers large fluctuations in the filopodia lengths as compared to in the absence of capping protein. 

The interplay between ratcheting and percolation has been explored to some extent in this work. Most percolation 
analysis has been reserved for "static" systems, where energetics dominates For example, recent Monte Carlo 
simulations of a system of binary rods, one rod to model the actin filaments, the other to model the cross- linkers, 
have explored the percolation transition in isotropic and bundled phases (53|. In the isotropic phase, the percolation 
transition resembles that of a vulcanization transition. It would be interesting to explore more thoroughly the 
percolation transition as a function of growing rods in the context of various phases ranging from the isotropic to 
dendritic to star phases as well as potentially new phases. 

Finally, star formation has recently been observed in noncquilibrium actin/cross- linker networks where polymeriza- 
tion is initiated [54[. These stars form in the absence of Arp2/3. We have seen this to be case in our calculations. 
Lieleg and collaborators [54| have measured the fractal dimension of the stars to be around 1.8. A study of the fractal 
dimension of our simulated stars has yet to be conducted. In any event, by keeping track explicitly of the G- actin 
monomers as the star forms, one enters the realm of diffusion limited aggregation (DLA) [55|. The angle-dependent 
branching and the bundling reorientation should lead to brand new cluster characteristics not yet observed in physical 
systems but, indeed, sought out by biological systems for reasons yet to be discovered. 
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Appendix A: Bundling due to bending 

To understand the role of bending, we construct the bundle size distribution for filaments emerging from an ordered 
two-dimensional array in which the filaments have the same initial angle of zero degrees with some regular spacing. 
Polymerizing filaments eventually bend in the transverse direction and "find" each other. It is the only way they 
can find each other. We compute the collision probability (or rate) for a semiflexible polymer to first reach some 
transverse distance, i.e. to run into another growing semiflexible polymer, denoted by Aj?™ d . As before, provided the 
angle is less then some critical angle, we assume realignment occurs. 

To compute A^ e ^ d , the Hamiltonian for a semiflexible polymer is given by 

where t(s) = dr ^ is the tangent vector of the curve r(s) and n is the bending rigidity, which is related to the 
persistence length, L p , in two dimensions via L p = [Ha] . Furthermore, the partition function is given by 

Z(t L ,L\t ,0)<x [ ( } L V[t(s)]S(\t(s)\-l)exp(—^-), (A2) 

Jt(0)=t k B± 



where the inextensibility condition is enforced by the delta function. Minimizing the semiflexible polymer paths yields 
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an equation for Z , namely, 

dZ 1 d 2 Z 



ds L p d0 2 



(A3) 



where 6(s) is the angle between t(s) and a fixed reference axis. In order to obtain information about the position 
vector, r, we replace the derivative on the left hand side of Eq. A3 with the convective derivative along the semiflexible 
polymer path, i.e. 

<Ts + ™®l + ^& z = ^W> (A4) 

Assuming 9 << 1 and integrating over x, since we arc only interested in the transverse fluctuations leading to 
bundling, simplifies Eq. A4 to 

Defining s' = s/L p and y' = y/L p yields a Fokkcr-Planck equation describing a particle undergoing random 
acceleration. In the absence of any boundary condition, this equation is easily solved. However, to compute a collision 
time one needs to consider this randomly accelerating particle in a finite interval < y < a where there are absorbing 
boundary conditions at either end of the interval. To date, only the steady state absorption/collsion probability has 
been calculated by Bicout and Burkhardt [57|. If we assume that 6q = and the filament to be bundled starts 
halfway in between two other growing filaments, then the steady state absorption probability is 1/2. Converting 
this absorption probability to a collision rates yields \\ e i d = 1 ; where p\(0) is the initial spacing between 

filaments and we have assumed a constant concentration of G-actin. Since we have averaged over x and retain only 
the y-coordinatc, we are considering a onc-dimensional concentration. To generalize A!j >c 1 nd to bundles of size i and k, 
L p (i) oc i for a loosely packed bundle such that 

A bcnd = A_ 1 Wl(0) + k -ly (Ag) 

hk 180° 2 L p y ' y ' 

This kernel is different from the case we analyzed in Section III because A b ™ d decreases with increasing i and k. The 
bundle pcrsistancc length increases as the bundle size increases making transverse fluctuations more difficult. 

For the above kernel, the existence of self-similar solutions for the kinetic aggregation equation, in accordance with 
the dynamical scaling hypothesis, has recently been proved [Hj]. In particular, 

Pi(t) ~ (-J=). (A7) 

Furthermore, f (x) is expected to scale as x~ a exp(Cx - ^) for small x such that Pi(t) is nonmonotonic in i, unlike 
kernel studied in Section III. We can test this result using the finite-dimensional simulations. 

In Sec. Ill we assumed previously that the cross-linkers are able to reorient, or bend, even the largest of bundles. 
However, there is presumably a cross-linker concentration-dependent bundle size beyond which the cross-linkers cannot 
reorient the bundle in spite of the small angle collision. As opposed to implementing a threshold bundle size for 
reorientation, one could implement a "soft" cutoff with A^fc ~ (Vi + -\/fc)(i _1 + fc _1 ). By incorporating this natural 
cutoff into the model, we can obtain a mean field bundle size distribution with a peak and an exponential tail. We do 
not have to propose a rate for filaments leaving the bundle, which is unlikely given that the filaments are extended 
objects, as was done in Ref. (26| . to obtain a peak. 
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FIG. 2: Log- linear plot of final bundle size distributions for different concentrations of Arp2/3 with 9 C = 20°. 
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FIG. 3: (a) The average bundle size as a function of the concentration of Arp2/3 with C — 20°. (b) The average bundle size 
as a function of C with a = 5nM. 




FIG. 4: (a) Scaling collapse for the intermediate bundle sizes with a = 10 nM, 8 C = 20°, and p m = 10 fiM . The curves clearly 
illustrate convergence towards the proposed scaling in this region, (b) Scaling collapse for large bundle sizes with a = lOnAf, 
0c = 20°, and p m = 10 fiM. 
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FIG. 5: (a) Snapshot of the simulation with Rt = 0.0 fim. The line in the lower righthand corner is 1 micron. The remaining 
parameters are: Nf = 500, Na = 250, 9 C = 30°, and Ri = 0.6, do = 0.003 \im, dx = 0.017 fJ,m, and dA = 0.01 /im. (b) 
Snapshot of the simulation with Rt, — 0.5 /im with the bead depicted in aqua (in color). The remaining parameters are the 
same as in (a) with the exception, Ri = 0.4 fim. 
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FIG. 6: Log-log plot the bundle size distribution at different times measured in simulation time steps. We use Nf 
N a = 250, R b = 0, R, = 0.6 pm, d A = 0.01 [im, d x = 0.017 /xm, Rent = 0.75 fjm, Q c = 30°, and d f = 0.03 fxm. 
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FIG. 7: (a) Log-linear plot of the bundle size distribution for different 8 c s with Nf = 500, Na = 250, Rb — 0, Ri = 0.6 (jm, 
dA — 0.01 fim, dx = 0.017 fim, R cu t — 0.75 \im, and df — 0.03 (jm. (b) Log-log plot of Fig. 7(a) The line has a slope of -0.5. 



20 




FIG. 8: Log-linear plot of the bundle size distribution for different Nas with Nf + Na = 750, Rt = 0, Ri = 0.6 (im, 
d,A = 0.01 fim, dx = 0.017 /im, R cut = 0.75 fj,m, 8 C = 30° and df — 0.03 /im. 
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FIG. 9: Q and the spatial distribution of the branched Arp2/3, Prob a , as a function of the radius for two different sets of 
parameters. One data set is with Nf = 100, Na — 650, Ri = 2 /ira and 8 C = 0. The second data set is with Nf = 100, 
N A = 650, Ri = 1.5 ^m, and 9 C = 0. 




r (microns) 



FIG. 10: Plot of Q(r) for the parameters plotted in Fig. [7] and the spatial distribution of branch points for 6 C — 0°, which is 
plotted as filled diamonds. 
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FIG. 11: Log-log plot the bundle size distribution at different simulation time steps. We use Nf = 500, Na = 250, Rt = 0.5 jim, 
i?i = 0.4 fj,m, d A = 0.01 nm, d x = 0.017 /jm, 7? cut = 0.75 /xm, 6> c = 30° and d/ = 0.03 fj,m. 
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FIG. 12: (a)Log-linear plot of the bundle size distribution for different 9 c s with Nf — 500, N a — 25 0, Rb = 0.5 /im, Ri = 1.2 /jm, 
(La = 0.01 /im, = 0.017 /xm, 7i cu t = 0.75 \im, and d/ = 0.03 fim. (b) Log-log plot of Fig. 12(a) The line has a slope of -1.0. 
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FIG. 13: (a) O(r) for different 8 c s with Rt = 0.5 )im and i?i = 0.4 /im. We also plot Prob a for the 



0° case, depicted in 



filled upside down triangles, (b) 0(r) for different (9 c s with Rt = 0.5 fj,m and Ri — 1.2 /^m and Prober for the # c = 0° case, 
again, depicted in filled upside down triangles. 
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FIG. 14: Snapshot of the simulation with Rt — 0.5 fj,m, Ri — 0.4 pm, 8 C = 30°, and c = 0.66 s 
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FIG. 15: Position of the filament tips as a function of the radius of the system for different capping rates. Here, Ri = 0.3 fj,m 
to accentuate the ratcheting effect. 




FIG. 16: Log-linear plot of the bundle size distribution for the ordered initial conditions. Here, Nf = 100, the filament seeds 
start 0.1 microns apart with the same initial angle, and dx = 0.017 fim. 



